…
library(lidR)# for working with LAZ files
library(sf) # for working spatial class
library(raster)# for working with raster
library(rayshader) # for 3D viz
library(rgl) # for interactive plots
motte_point = st_read("data/canmore_point.shp", quiet = TRUE)
motte_buffer = st_buffer(motte_point,dist = 50)
Next, I am reading the LAZ files and clipping the point cloud to the extent of immidiate area around the motte.
#read laz files
las = readLAS("data/NX6055_4PPM_LAS_PHASE3.laz")
motte = clip_roi(las, motte_buffer)
#plot lidar point cloud
plot(motte, bg = "white")
rglwidget()
grab and rotate the model !
In this step I am running standard algorithms from LiDR package to compute Digital Surface Model (DSM) and Digitial Terrain Model (DTM).
rgl.clear()
# create dsm and dtm
dsm = grid_canopy(motte, 0.5, pitfree())
# assign coordinate system
crs(dsm) = CRS("+init=epsg:27700")
# create dtm
dtm = grid_terrain(motte, 0.5, algorithm = knnidw(k = 6L, p = 2))
# addign coordinate system
crs(dtm) = CRS("+init=epsg:27700")
par(mfrow = c(1,2))
plot(dtm, main = "DTM", col = height.colors(50))
plot(dsm, main = "DSM",col = height.colors(50))
# dsm
slope_dsm = terrain(dsm, opt = 'slope')
aspect_dsm = terrain(dsm, opt = 'aspect')
hill_dsm = hillShade(slope_dsm, aspect_dsm, angle = 40, direction = 270)
# dtm
slope_dtm = terrain(dtm, opt = 'slope')
aspect_dtm = terrain(dtm, opt = 'aspect')
hill_dtm = hillShade(slope_dtm, aspect_dtm, angle = 5, direction = 270)
#plot
par(mfrow = c(1,2))
plot(hill_dtm, main = "DTM Hilshade", col = grey.colors(100, start = 0, end = 1),
legend = FALSE)
plot(hill_dsm, main = "DSM Hillshade", col = grey.colors(100, start = 0, end = 1))
#And convert it to a matrix:
elmat = raster_to_matrix(dtm)
elmat %>%
sphere_shade(texture = "imhof1") %>%
add_shadow(ray_shade(elmat, zscale = 0.5, sunaltitude = 30,lambert = TRUE),
max_darken = 1) %>%
add_shadow(lamb_shade(elmat,zscale = 0.5,sunaltitude = 30), max_darken = 0.2) %>%
add_shadow(ambient_shade(elmat, zscale = 0.5), max_darken = 0.2) %>%
plot_3d(elmat, zscale = 0.5,windowsize = c(1000,1000),
phi = 40, theta = 180, zoom = 0.8, fov = 1)
rglwidget()
grab and rotate the model !